{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Further Hypothesis Testing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Select this cell and type Ctrl-Enter to execute the code below.\n",
    "\n",
    "library(tidyverse)\n",
    "\n",
    "set_plot_dimensions <- function(width_choice, height_choice) {\n",
    "    options(repr.plot.width=width_choice, repr.plot.height=height_choice)\n",
    "}\n",
    "\n",
    "cbPal <- c(\"#E69F00\", \"#56B4E9\", \"#009E73\", \"#F0E442\", \"#CC79A7\", \"#0072B2\", \"#D55E00\")\n",
    "\n",
    "set_plot_dimensions(5, 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# You should see \"Attaching packages\" and some ticks by the packages loaded.\n",
    "# The \"Conflicts\" aren't a problem.\n",
    "\n",
    "# Other problems loading the library? Try running this cell.\n",
    "\n",
    "install.packages('tidyverse')\n",
    "\n",
    "library(tidyverse)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7 - Correcting for multiple hypothesis tests"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run this cell to load the data.\n",
    "\n",
    "data <- read_csv(\"stars.csv\")\n",
    "\n",
    "type_key <- c('Brown Dwarf', 'Red Dwarf', 'White Dwarf', 'Main Sequence', 'Supergiant','Hypergiant')\n",
    "spectral_classes <- c('O','B','A','F','G','K','M')\n",
    "\n",
    "data$type <- factor(data$type)\n",
    "data$spectral_class <- factor(data$spectral_class, levels=spectral_classes)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unfortunately, there is a problem with the previous analysis.\n",
    "\n",
    "Recall that the significance level, $\\alpha$, is defined as the probability of incorrectly rejecting $H_0$ when it is actually true (i.e. the probability of a Type I error).\n",
    "\n",
    "When we perform [*multiple related hypothesis tests*](https://en.wikipedia.org/wiki/Multiple_comparisons_problem), we increase the chances of producing such a Type I error.\n",
    "\n",
    "For example, if $\\alpha=0.05$ and we perform 100 tests, we *expect* to generate 5 Type I errors. This can be a serious problem when large numbers of hypothesis tests are carried out simultaneously, for example in screening thousands of genes for association with a disease.\n",
    "\n",
    "We therefore need a strategy to control the rate of Type I errors. A very simple approach is given by the [*Bonferroni correction*](https://en.wikipedia.org/wiki/Bonferroni_correction):"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bonferroni correction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Theory\n",
    "\n",
    "When conducting $n$ related hypothesis tests, we reduce the significance level for each test to $\\alpha/n$.\n",
    "\n",
    "The probability of making a Type I error *over the whole set of tests* (known as the *family-wise error rate*, FWER) therefore remains at $\\alpha$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Application"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Shapiro-Wilk test for normality\")\n",
    "print(\"\")\n",
    "\n",
    "p_values <- c()\n",
    "alpha <- 0.05\n",
    "n <- 6\n",
    "\n",
    "for(t in seq(0,n-1)){\n",
    "    sample <- data %>%\n",
    "        filter(type == t) %>%\n",
    "        pull(temperature) %>%\n",
    "        log\n",
    "    p_values <- append(p_values, shapiro.test(sample)$p.value)\n",
    "}\n",
    "\n",
    "print(paste(\"with uncorrected alpha =\",signif(alpha,3),\":\"))\n",
    "for(i in seq(1,n)){\n",
    "    result <- \"\"\n",
    "    if(p_values[i] < alpha) result <- \"*** REJECT H0 ***\"\n",
    "    print(paste(i, type_key[i], \": p =\", signif(p_values[i],3), result))\n",
    "}\n",
    "\n",
    "print(\"\")    \n",
    "\n",
    "print(paste(\"with Bonferroni correction, alpha/n =\",signif(alpha/n,3),\":\"))\n",
    "for(i in seq(1,n)){\n",
    "    result <- \"\"\n",
    "    if(p_values[i] < alpha/n) result <- \"*** REJECT H0 ***\"\n",
    "    print(paste(i, type_key[i], \": p =\", signif(p_values[i],3), result))\n",
    "}\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "After correcting for multiple hypothesis testing, the red dwarf p-value is not significant.\n",
    "\n",
    "We should report to Professor Xu that log(temperature) is not normally distributed for the brown dwarf, supergiant and hypergiant types.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Alternative methods for multiple testing correction\n",
    "\n",
    "The Bonferroni correction is simple to apply, but it may be too conservative when there is a very large numbers of tests, or when the tests are not independent (for example, genes are often related to other genes so are likely to share properties).\n",
    "\n",
    "The [*Benjamini-Hochberg procedure*](https://en.wikipedia.org/wiki/False_discovery_rate#Benjamini–Hochberg_procedure) is an alternative approach. Instead of controlling the FWER, this method controls the *proportion of the positive tests that are incorrect*, i.e. the proportion of rejected $H_0$'s that are Type I errors. This is known as the *false-discovery rate*, FDR."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If a vector of p-values is available, the `p.adjust()` function will compute the *adjusted* p-values according to the Benjamini-Hochberg method (or several other available methods). Adjusted p-values are sometimes called *q-values*. These can then be compared to the original $\\alpha$ value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# In our example, the Benjamini-Hochberg method is less conservative than Bonferroni: \n",
    "# The red dwarf p-value still appears to be significant when we use this method.\n",
    "\n",
    "q_values <- p.adjust(p_values, method='BH')\n",
    "\n",
    "print(paste(\"with Benjamini-Hochberg correction, alpha =\",signif(alpha,3),\":\"))\n",
    "for(i in seq(1,n)){\n",
    "    result <- \"\"\n",
    "    if(q_values[i] < alpha) result <- \"*** REJECT H0 ***\"\n",
    "    print(paste(i, type_key[i], \": q =\", signif(q_values[i],3), result))\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
